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ABSTRACT 



>■ 

T-H Aims. We numerically explore electron acceleration and coronal heating by dissipative electric fields. 

^^ Methods. Electrons are traced in linear force-free magnetic fields extrapolated from SOHO/MDI magnetograms, endowed with anomalous 
^^ resistivity (;;) in localized dissipation regions where the magnetic twist V x b exceeds a given threshold. Associated with ;; > is a parallel 
>D electric field E = z/j which can accelerate runaway electrons. In order to gain observational predictions we inject electrons inside the dissipation 
regions and follow them for several seconds in real time. 

Results. Precipitating electrons which leave the simulation system at height z = are associated with hard X rays, and electrons which escape 
O i' at height z ~ 3- 10* km are associated with normal-drifting type Ills at the local plasma frequency. A third, trapped, population is related 
to gyrosynchrotron emission. Time profiles and spectra of all three emissions are calculated, and their dependence on the geometric model 
parameters and on 77 is explored. It is found that precipitation generally preceeds escape by fractions of a second, and that the electrons perform 
many visits to the dissipation regions before leaving the simulation system. The electrons impacting z = reach higher energies than the 



o 

Q 
I 

O 

jrt , escaping ones, and non-Maxwellian tails are observed at energies above the largest potential drop across a single dissipation region. Impact 

• • • maps at z = show a tendency of the electrons to arrive at the borders of sunspots of one polarity. 
. 5^ ■ Conclusions. Although the magnetograms used here belong to non-flaring times, so that the simulations refer to nanoflares and 'quiescent' 
K^ , coronal heating, it is conjectured that the same process, on a larger scale, is responsible for solar flares. 

5t 1 '^^y words. Acceleration of particles - Sun: corona, flares. X-rays, radio radiation - Methods: numerical 

1 . Introduction cases where the type III onset could not properly be resolved. 

At millimeter and decimeter wavelengths, synchrotron emis- 

Observations of radio waves and hard X-rays (HXR) from . , r^- - ,1 i ^rvrv^ t j • . • j- 

■' sion (e.g., Gimenez et al. I2005L decimetnc radio continuum 

solar flares allow the study of the acceleration and propaga- j j • . • 1 .■ ^c ■ » tt-i • o r> Hnm t j 

■' ±- ±- o ^jj(j (jecimetric pulsations (Saint-Hilaire & Benz I2003L and 

tion of high-energy electrons which are responsible for both , . ^ . ., //^ •• j 1 . 1 haah i, i, 1 u c j .. 

° °-' ^ decimetnc spikes (Gudel et al. 119911 1 have also been found to 

types of emission. The picture emerging from the observa- , • . j -^i, uvr. 1,1 • .1, u j 

■'^ r & b ^Q associated with HXR, roughly in the above decreasing or- 

tions is generally complicated (Benz et al.y005). While HXR j ^ • .• u ui-. -ru . ^ j • .• j 

° ■' ^ ^ ' ' der of association probability. These types of decimetnc and 

(bremsstrahlung) emission IS often associated with metric type .„. .. • j- • ■ 1 . j» »• n » j 

^ °^ ■'^ millimetric radio emission are related to magnetically trapped 

III radio groups (Aschwanden et al. .19953 , the timing of in- , , i, 1 • . ui-.- u i- j » u 

° ^ ^ ' ' ° electrons, where loss-cone instabilities are believed to be re- 
dividual type Ill's and HXR fine structure is erratic. Metric ui ^ . 1 . . /r^ ■■ o ci *.- i im/; i 
-'^ sponsible tor temporal structures (Kuijpers & Slottje 119761 

type Ill's are presumably caused by electron beams exciting Aschwanden & Benz[Mi Fleishman & Melnikov[MH». 
Langmuir waves, which then couple to electromagnetic (ob- 

, , N ■ e t- »u • f , . u The variety of observed behaviour accounts, on the one 

servable) modes. Sometimes, there is perfect agreement be- ^ 

jjj . , TjvTj • u ^- ^u *u hand, for the non-linear (coherent) radio emission processes, 

tween type III onsets and HXR maxima, but in other cases there ^ ' ^ 

, . , , , , ,. n^, . . On the other hand, it also reflects the geometric complexity 

IS no obvious peak-to-peak correlation, or one of the emissions ° f j 

■ , . ,, ., . . A tu t TTT' t e A of the active regions, which results from the nonlinear dynam- 

is absent altogether. As a trend, the type III s onsets were found ° -' 

»uji juf*- c A tti, TTvr. c ics in the convection zone. In fact, the complex behaviour and 

to be delayed by fractions of a second against the HXR fine ^ 

structures (Aschwanden et A^ Arzner & Benz^OQSaj, al- *^ '^""^P^^'^ 8^°"^^^^ ^'^ ^^^^^ t° ^^ '^^"^^^ ^"'^ ^^^""'^ °f 

., 1, 4.1, 4, TTT j^ A t^ • . I. A ^ ■ each Other. Including a realistic amount of geometric complex- 
though the type III frequency drift may assist such delays in ° fe f 
ity into a numerical model of solar flares was a major motiva- 

Send ojfprint requests to: K. Arzner tion for the present work. 
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Apart of geometrical aspects, the modeling of solar flares 
requires the specification of a physical acceleration mecha- 
nism. We shall assume here, as a working hypothesis, that ac- 
celeration is caused by DC parallel electric fields due to anoma- 
lous resistivity. Such fields are capable to accelerate high- 
energy electrons for which the electric force is no longer coun- 
terbalanced by the collisional drag (Dreicer 1960 1. It should be 
pointed out that runaway acceleration applies only to a small 
fraction of electrons in the high-energy tail of the electron dis- 
tribution, and that the majority of electrons is (Joule) heated 
rather than accelerated to superthermal energies. The electric 
fields envisaged above are of macroscopic and dissipative na- 
ture, and mark the irreversible release of non-potential mag- 
netic energy. We thus follow a scenario, proposed by Parker 
<1972ll983ll993l l. where the random photospheric footpoint 
motion twists and shears the coronal magnetic loops such as 
to develop tangential discontinuities. This happens intermit- 
tently and in multiple localized regions throughout the solar 
corona, and does not generally imply loop instability or global 
re- structuring of the solar corona. Parker's idea has initiated 
many investigations on coronal heating and -eruptions (e.g., 
Schumacher & Kliem 1996 Torok & Kliem 2001 Gudiksen 
& Nordlundl2002i, and has significantly improved the under- 
standing of the global flaring process. In our view, a strong but 
heuristic argument for adopting dissipative (rather than con- 
servative) electric fields as responsible for particle acceleration 
in solar flares is the intermittent and violent nature of flares, 
hinting to a catastrophic process which could not be cast in a 
Hamiltonian formalism. However, conservative electric fields, 
as used in most investigations on stochastic acceleration (e.g. 
Karimabadi et al. 1987, Schlickeiser ■2003J . may also act as 
particle accelerators. 

In the present study, we adopt a test particle approach 
similar as in the simulations of Matthaeus & Lamkin ( 1986 1, 
Dmitruk et al. (2003 2004 1, Arzner & Vlahos (2004i, and 
Arzner et al. (2006 1, but with two important modifications. 
First, the electromagnetic force fields are not taken from MHD 
turbulence simulations or random-phase turbulence proxies, 
but from observed magnetograms. Secondly, in order to span 
the many orders of magnitude between the electron Larmor ra- 
dius and the size of magnetic loops, we allow for gyrokinetic 
motion, keeping track of the gyro phase in an approximate way. 
This technique enables us to follow the electrons over several 
seconds in real time (several 10** gyro times), thus reaching the 
time scales on which solar flares are observed. In the simula- 
tions, runaway electrons are injected at f = inside the local- 
ized dissipation regions, and followed numerically along (and 
sometimes across) the magnetic field lines. Observational pre- 
dictions for HXR and radio waves are obtained as the electrons 
impinge the chromosphere, get trapped, or escape to the higher 



The article is organized as follows. Section 2 describes the 
construction of the coronal field; Section 3, the particle orbits; 
Section 4, the numerical results and observational predictions. 
These are then summarized and discussed in Section 5. 



2. Coronal electromagnetic fields 

Our domain is a slab of height //, bounded at z = by the 
photosphere, and filled with time-independent force-free (e.g., 
Gary '1989^ magnetic fields. We do not ask here for a perfect 
field reconstruction, but for a generic configuration which is 
compatible with the observation. Accordingly, we use an ele- 
mentary linear force-free extrapolation from the normal photo- 
spheric magnetic field. This represents a local approximation 
at best (Wheatland I1999> . and we have a slab thickness H of 
not more than a few 10"^ km in mind when fixing the physical 
scaling. (The linear force-free assumption predicts loops flar- 
ing with height, in contrast to the slender high-ranging loops 
observed by TRACE.) Further characteristic length and time 
scales are compiled in Table [0 



Symbol 


Meaning 


Typical value 


n-' 


particle time unit 


10-^ s 


/o 


particle length unit 


0.3 m 


/ 


magnetogram resolution 


660 km 


L 


magnetogram size 


l.T-lO^km 


H 


slab height 


3 • 10* km 


I /a 


force-free scale 


2 ■ 10* km 


l/"c 


critical twist scale 


3 ■ 10^ km 



Table 1. Notation and physical scaling. 



2.1. Construction of the force-free magnetic field 

The force-free condition requires that V x B = aB, where a 
is assumed to be constant ('linear' force-free field). Our con- 
struction of the force-free magnetic field follows the lines of 
Alissandrakis ( 1981 1, with modifications concerning the selec- 
tion and number of Fourier modes. In order to avoid large but 
passive 3-dimensional arrays, the strategy is to locally com- 
pute B(x) from a restricted number of Fourier modes b(k) e'*'". 
In terms of these, the force-free condition becomes /k x b(k) - 
Q'b(k). This equation has a non-trivial solution only if |kp - oP-, 
which we write in the form 



/t. = + ^a^ - k]_ where k^ = kl - k^ ., 



and if b(k) is proportional to 

(—k^k, + iaky, —kyk^ — iak^, k^ + k ) . 



(1) 



(2) 



See MacLeod (I1995> for a deeper discussion of the curl 
eigenfunctions. We assume kx and ky to be real, so that the 
magnetic field in the (x, y)-plane is bounded, while kr may be 
complex with positive imaginary part. Thus modes decaying or 
oscillating with height are permitted, and we shall use a small 
admixture of oscillating modes to transport structures from the 
chromosphere to larger heights. From a physics point of view 
this procedure is justified by noticing that the outstreaming so- 
lar wind (in a Id slab geometry) does not require B(z) to vanish 
as z — > oo. Furthermore, diff'erent wave vectors are not ratio- 
nal multiples of each other, so that the magnetic field extends 
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normalizing the eigenvectors such that \b^\ - 1, and using 
the (approximate) orthogonaHty of the harmonic functions over 
Lx L, the coefficients C(k) are given by 



C(k) 



20 40 60 80 100 120 
x/l 

B.(x,y,0) 



50 -SO '50 EnC 250 




20 40 60 80 100 120 

x/l 



50 100 150 2C(I 250 



Fig. 1. Boundary fields (top) and their sparse-Fourier rep- 
resentations (bottom) of a bipolar configuration (left) and 
a SOHO/MDI magnetogram (right) recorded on August 17, 
2002. 10:40 UT. The lower left corner of the magnetogram 
corresponds to (78", -36") in heliocentric coordinates, and the 
scale is / = 0.94". 





Fig. 2. Power spectral density (grayscale) of the boundary fields 
of Fig. [T] (top). Left: bipolar configuration; right: SOHO/MDI 
magnetogram. The sparse-Fourier components are marked by 
dots. 

in a non-periodic, generic, way across the (x,y)-plane. The ir- 
regular spacing of the Fourier components also avoids aliasing 
artifacts, and their small number allows a non-expensive eval- 
uation of the magnetic field even if no fast Fourier transform is 
available. However, the restriction to few Fourier components 
makes an exact matching of the boundary conditions impossi- 
ble. Instead, we require agreement in a least square sense only, 
and minimize the mean-square deviation of the sparse-Fourier 
field B^(x, y, 0) and the photospheric boundary field B,o(x, y) in 
the square Lx L. Setting 



B(x) = Re ^ Ckb(k)e''"' , 



const. 



/ 



dxdy fo*(k)e 

LxL 



—ik^x—ik^y 



Bzo(x,y). 



(4) 



The boundary field B-o(x,y) is given on a square lat- 
tice (cell size f), and the integral @ is computed numeri- 
cally. 'Complementary' solutions of the homogeneous bound- 



(3) 



ary value problem (Chiu & Hilton 119771 Petrie & Lothian 
|2p03) are disregarded. The sparse-Fourier wave vectors (k^, ky) 
are selected by the following procedure. First, the fast Fourier 
transform of the discrete boundary field is computed and those 
Nii (regularly spaced) wave vectors are chosen which have 
largest power spectral density. Then, a random perturbation 
{Ak^, Aky) :$ n/l is added in order to break the exact periodicity. 
The eigenvectors (Eq.|3 automatically ensure that V ■ B = 0. 

A test case with two Gaussian footpoints of opposite polar- 
ity is shown in Fig. Q] (left column), with the true boundary 
field B,Q(x,y) presented in the top left panel and its sparse- 
Fourier version B,(x, y, 0) presented in the bottom left panel. 
The sparse-Fourier version contains Nk = 100 components 
only. This low number already provides a good approxima- 
tion (correlation coefficient 0.94), because the Fourier spectrum 
(Fig. 121 left) is concentrated at the origin. When real data are 
used, a larger number of Fourier components is needed. This is 
illustrated in Figure ^ (right column), showing a SOHO/MDI 
magnetogram (top right) and its sparse-Fourier representation 
(bottom right) using 5000 Fourier components (Fig. |2] right). 
The magnetogram is located close to the center of the so- 
lar disc, so that projection effects are negligible. The mag- 
netogram resolution is / - 660 km (rescaled from the origi- 
nal SOHO/MDI resolution), and the correlation coefficient be- 
tween the magnetogram and its sparse-Fourier representation is 
0.93. The magnetogram belongs to a non-flaring configuration 
and therefore the force-free extrapolation should be a reason- 
able approximation. Although the smallest scales are not well 
represented by the sparse-Fourier approximation, they rapidly 
decay with height (Eq.[T), so that the approximation becomes 
better as z increases. This fact is used to save computation 
time by restricting the field computation to components with 

U'V'I > 10-3. 



2.2. Properties of the magnetic field 

Figure|3]displays the force-free magnetic fieldlines of the bipo- 
lar test case (top) and the SOHO/MDI magnetogram (bottom). 
The a parameter is such that al = 0.03, corresponding to a 
force-free scale 1/a of about 2-10'* km. The field lines start out 
at z = at random with density proportional to |B|. Several field 
line integration schemes have been tested, and the force-free 
condition was also verified using finite 3D differences. Most 
field lines connect the two poles of opposite polarity, but not 
all, because the sparse-Fourier field does not vanish outside the 
poles, nor outside LxL. 

The absolute scaling of B is constrained by the SOHO/MDI 
data which give, for Figure |3] (bottom), a longitudinal pho- 
tospheric rms field -\J{Bzo(x,y)^) - 50G, with excursions to 



K. Arzner and L. Vlahos: Gyrokinetic electron acceleration 





^^Q 



Fig. 3. Sparse-Fourier force-free magnetic field lines extrapo- 
lated from bipolar (top) and magnetogram (bottom) boundary 
data. 



+400G at the footpoints. These values are typical for active 
regions; the average over the quiet solar surface is smaller 
(~1G, see Kotov "2002). The mean-square magnetic field 
strength is related to the Fourier amplitudes by {B^(x, y, 0)^) = 
4 2k I^A(k)P and similar relations for B,, and B^. 



2.3. Electric field 

Since the configuration is assumed to be time-independent, 
the electric field has only an Ohmic component E = 77J. 
The (scalar) resistivity 77 is either zero or anomalous, i.e., 
much larger than the collisional value T/spitzei (Helander 2002 1. 
Regions with 77 » T/spitzei are called here 'dissipation regions'. 
Owing to the force-free condition V x B = aB, the electric field 
E = rjaB is purely parallel. 

Different criteria for the occurrence of anomalous resistiv- 
ity have been proposed, relying on linear micro-instabilities 
or heuristic arguments. A well-known argument limits the 
electric current to jc - encs where c, is the speed of sound 




1^^ 



T,'o'^ 



Fig. 4. Different definitions of dissipation regions. Top: current 
threshold. Bottom: magnetic twist threshold. Both panels refer 
to the magnetogram of Fig.|3](bottom). It is the twist threshold 
which is used in the sequel of this paper. 

produce shocks and dissipation in the cuiTent-caiTying electron 
fluid. The resulting dissipation region |j| > jc is depicted in Fig. 
0](top) for the configuration of Fig. |3l (bottom). As the force- 
free cuiTent is proportional to the magnetic field, the dissipation 
regions simply delineate the magnetic field strength. 

Another approach invokes the twist or shear of magnetic 
field lines rather than the electric current density. The physical 
motivation for this stems from the observation that regions of 
larger magnetic field should also be able to guide larger cur- 
rents before disruption. Thus the automatic increase of j^ with 
increasing |B| should be compensated, which is achieved by 
considering the quantity 



^ = 1 



(5) 



rather than B. The criterion for the occurrence of anoma- 
lous resistivity is then defined by 



IV X bl < Ur 



(6) 



(Papadopoulos [T980l l: an excess |V x B| > jc is then argued to 



The resulting dissipation regions are shown in Fig. |3 (bot- 
tom); the top and bottom panels represent equal dissipation vol- 
umes. 

It foUows from the force-free condition that b ■ V x b = a. 
Thus, for constant a, the field-aligned component of the mag- 
netic twist V X b is constant throughout the whole volume. This 
implies, in particular, that |V x b| can never be smaller than 
a, and that the threshold Uc must be larger than a. In Fig. |4] 
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(bottom), Ucia - 6.6. Since small scales with k^ > a de- 
cay exponentially with height, the steep gradients, and thus 
the dissipation regions, tend to accumulate at low altitudes. 
This is in agreement with the observation that the coronal 
heating function is localized within a height range :$ 10"^ km 
( Aschwanden '200 1 V The dissipation regions are thus not uni- 
formly distributed in height, and the dissipative volume Vd can 
be characterized by a column height 



hd = VdlL^ 



(7) 



rather than by a volume fraction referring to the full numer- 
ical domain L^H. In the case of Fig.l^one has hd - 0.19/. 

The condition (|6jl is local. Non-local versions of Eq. (|6jl 
have been used by Georgoulis & Vlahos ( 1998) along pairs of 
magnetic field lines, and by Vlahos & Georgoulis ( 2004 1 on a 
discrete lattice. We use here Eq. (|6|l, since it can be evaluated 
at the particle position and thus simplifies the simulation. The 
threshold Uc is chosen such that uj - 0.2. 

3. Particle orbits 

3.1. Dimensionless units and physical scaling 

The system of units used in this article is similar as in Arzner & 
Vlahos J2004> . Time is measured in units of Q ' where Q^ = 
5 2k |b(k)p is the non-relativistic mean-square gyro frequency 
at z = 0. Velocity is measured in units of the speed of light; 
energy is measured in units of mc^; momentum is measured in 
units of mc; distance is measured in units of Iq = cQ.\ We 
take <|B(x,y,0)|2)'/2 = 50G (see Sect.E3 to fix the absolute 
scaling, so that the electron cyclotron frequency is 140 MHz, 
the time unit is Q"' - 10"''s, and the length unit is Iq - 0.3m. 
The electric field is measured in units of cB, and the resistivity 
77 is measured in units of c^Q ' . 

3.2. Exact orbits 

In the above unit system, the exact particle equations of motion 
read 



dx 
dt 

y— = V X B + E - (v ■ E)v 
dt 



(8) 
(9) 



where E = z/V x B and y - 1/ Vl - v^ is the Lorentz fac- 
tor. Equations ^ and (|9jl are integrated by a traditional Cash- 
Karp/Runge-Kutta method (Press et al. I1998> with adaptive 
time step. The sparse-Fourier fields are computed at a lower 
rate, together with their Jacobians, from which a linear extrap- 
olation to the actual particle position is made. The time step is 
chosen such that both the Cash-Karp error and the field extrap- 
olation error remain within given bounds (AppendixIXt. 

3.3. Gyrokinetic approximation 

In the magnetic fields considered here, ranging from a few to 
a few hundred Gauss, and for typical kinetic energies between 



1 keV and 1 MeV, the electron Larmor radius -yv^/Q. varies 
from centimeters to hundreds of meters, and is therefore small 
compared to the resolution of the magnetic field (/ ~ lO"' km). 
The motion is thus almost always adiabatic, with possible ex- 
ceptions at critical points of the magnetic field, and inside the 
dissipation regions. This suggests a gyrokinetic approximation, 
with the possibility to switch to exact orbits if necessary. 

We base our guiding center mover on a relativistic (Brizard 
[1999") version of the gyrokinetic Equations of Littlejohn (il981l 
^983). The variables considered are (X, /:i||,0,/i), where X is 
the guiding center position, p^ - yvy is the momentum parallel 
to the magnetic field, is the gyro phase, and /i = j^y^v^/B is 
the magnetic moment, which is a motion integral in the present 
approximation. The coordinates (X,/?||) evolve according to 



dt B*\y y I 

^.-lB-.(^vB-rl 

dt B* \y I 

where y - Jl + p^^ + 2jiB and 



(10) 
(11) 



B* = B+p^Vxb 
db 



E* 



E-Pil 



dt 



(12) 
(13) 



with B* = B* ■ b. For our time-independent force-free mag- 
netic field, E* equals E and the last term in Eq. jlOt vanishes. 

In order to switch back to exact orbit integration we need 
to keep track of the gyro phase. In gyrokinetic approximations, 
the gyro phase is not unambiguously defined but subject to 
gauge freedom. Accordingly, various definitions have been pro- 



posed. We follow here Littlejohn J1988> . Let e be a unit vector 
in direction of the initial (f=0) Larmor radius. The vector Ci is 
perpendicular to the direction b of the local magnetic field, and 
a local cartesian triad is completed by setting §2 = Ci x b. The 
task is to follow the unit vector Ci along the gyro center, in- 
troducing a minimum of twist. This is achieved by solving the 
parallel transport (Fermi- Walker type) Equation 



dCi 
~dl 



-<t-) 



(14) 



where s is the distance along the gyrocenter orbit. The gyro 
phase, with respect to the direction Ci, is then given by = 
L y"'|B(X(f'))l<^f', and the exact orbit parameters (x, v) can 
be retrieved from (X, v\\,p, 0, Ci). More details and the bench- 
marking of our numerical approach are described in Appendix 

El 

4. Simulation results and observational 
predictions 

The slab-shaped simulation domain suggests a simple scheme 
of observational diagnostics, which is sketched in Figure [S] In 
this scheme, HXRs are associated with electrons impacting z - 
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0, where they are stopped by the rapidly increasing density ^ 
Electrons which remain magnetically trapped are associated 
with gyrosynchrotron emission, and electrons which leave the 
simulation domain at z = // are associated with radio type III 
emission at the local plasma frequency. 



" RADIO TYPE III 




z = H 



z = 



HXR 



Fig. 5. Schematic allocation of different emissions from high- 
energy electrons. Precipitating electrons are associated with 
HXRs; escaping and trapped electrons are associated with radio 
emission. Black lines symbolize magnetic field lines; boldface 
ellipses symboUze dissipation regions. 

We have performed several simulation runs with varying 
parameters of the electromagnetic fields. One of these runs is 
discussed in detail in Sections 14.11 to 14.41 while the outcome 
of the others is summarized in Section |431 In all simulations, 
electrons with initial velocity vo = 0.1c (£0 = 2.6 keV) are 
injected at f = inside the dissipation regions, and traced until 
one of the following conditions is met: z < ('precipitation'), 
z> H ('escape'), or f > T^ax ('trapping'). The simulation du- 
ration Tmax is long compared to the average arrival times at 
z = and z-H. r^ax is also long compared to the free escape 
time T() = HJVQ calculated on grounds of the initial velocity. 
The synchronized injection at f = is somewhat artificial since 
we use magnetograms of a non-flaring active region. We have 
chosen this initial condition in order to model HXR and radio 
transients and thus provide additional (timing) diagnostics in 
situations where an observable flare pulse occurs, based on the 
assumption that flares are large-scale versions of the process 
considered here. The predictions will, though, be qualitative 
only, since the time evolution of the flaring region is not taken 
into account. 

4.1. Electrons 

The outcome of a typical simulation is summarized in Fig.|6l 
The top panel represents the exit rates at z = and z-H. Most 
particles exit through the z = boundary and the precipitation 
rate has a sharp peak at Qf ~ 10^ (not resolved in Fig. |6j. 
The escape rate peaks later. The remaining (trapped) popula- 
tion gradually decreases (Fig.|6lmiddle; times smaller than the 

' Strictly speaking, z = delineates the photosphere probed by the 
SOHO/MDI magnetogram, and not the density step (transition region) 
seen by the precipitating electrons. We shall not make this distinction 
here since our model does not include a background density profile. 




Fig. 6. Precipitating, escaping, and trapped electrons. Top: exit 
rates at the slab boundaries. Middle: remaining (trapped) par- 
ticles. Bottom: terminal kinetic energies. Q.t - 5 ■ 10^ corre- 
sponds to 5 seconds in real time. 
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Fig. 7. Simulated HXR and type III light curves. 



free escape time Qtq = 6 ■ 10^ are not shown). A detailed look 
on the particle orbits shows that these usually perform numer- 
ous visits to the dissipation regions before reaching relativis- 
tic energies. The fraction of time spent by the electrons inside 
the dissipation regions is quite different for the different popu- 
lations. While precipitating particles spend about 10% of their 
time inside the dissipation regions, the escaping ones do so only 
for about 2%. The trapped ones spend 14% of their time inside 
the dissipation regions, but do not systematically gain energy. 
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In contrast, the escaping and - more pronounced - precipitating 
particles gain relativistic energies. The terminal energy spectra 
are shown in Fig. |6] (bottom). Non-Maxwellian tails occur in 
the precipitating spectrum at energies Ejmc^ ~ 1 . The trapped 
electron spectrum (Fig. |6l light gray) refers to the end of the 
simulation (t - Tmax)- 
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Fig. 8. Thick-target bremsstrahlung spectrum from the precipi- 
tating electrons. 



4.2. Hard X rays 

Light Curve and Spectrum. In our model, the precipitating 
population is associated with HXR bremsstrahlung, with the 
HXR time profile being directly proportional to the electron 
exit time density at z = (Fig. 0. Since z = is an ab- 
sorbing boundary we use the thick-target approximation. The 
(orientation-averaged) bremsstrahlung cross section is taken 
from the fully relativistic formula 3BN of Koch & Motz (1959 1. 
The resulting HXR spectrum is shown in Figure|8] At energies 
up to ~40 keV, it can be fitted by a (hard) power law with in- 
dex 1.4 (Fig. |8] inlet). At higher energies, it decays exponen- 
tially and finally drops oflF faster than exponentially at ~800 
keV. Most of the exponential part and the super-exponential 
decay could usually not be observed with X-ray observatories 
like RHESSI because of the limited dynamic range (2-3 orders 
of magnitude for M class flares; see Grigis & Benz 120041) . 

Impact Map. The present simulation also predicts the sites at 
which electrons impinge onto the chromosphere (Fig. |9] left), 
and allows a comparison of this with the boundary magnetic 
field (Fig. fright, repeated here from Fig.^for better clarity). 
Since we inject particles in all dissipation regions simultane- 
ously, the impact map represents a statistical prediction. There 
are two observations to be made: first, the electrons preferen- 
tially precipitate at negative magnetic polarity. Secondly, the 
impact density peaks not at the center of the sunspots or fil- 
aments, but at their borders. The first observation is a simple 
consequence of the constant-a assumption, by which the elec- 
tric acceleration is always parallel or antiparallel to the mag- 
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Fig. 9. Electron precipitation sites (left) and boundary magnetic 
field (right). Positive magnetic polarity is light; negative polar- 
ity is dark. 
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Fig. 10. Gyrosynchrotron intensity according to the exact ex- 
pression (Eq. [21 - crosses) and its approximation (Eq. [^ - 
dotted line). 

netic field^. The second observation is a consequence of the 
mirroring in converging magnetic field lines, which makes it 
hard for electrons to penetrate down to z = in regions of max- 
imal |B(ji:,3;, 0)|. The actual impact map represents a trade-off 
between the number of downstreaming electrons and their re- 
flection probability, which is optimal at the sunspot/filament 
boundaries. 

4.3. Electron beams 

Electrons which escape to altitudes z > H are considered as 
a proxy for type III radio bursts. Contrary to the incoher- 
ent bremsstrahlung and gyrosynchrotron radiation, the plasma 
emission of type III bursts is coherent and therefore not pro- 
portional to the number of emitting electrons. The non-linear 
dependence (say, oc N^ with ^ > 1) compensates for the rela- 
tively tenuous escaping population. The resulting time signal 
is shown in Fig. together with the HXR light curve. The 
scaling of the two intensities is arbitrary. As the background 
plasma density is not specified in our model, we cannot pre- 



- We have assumed here positive a and that the electrons have pos- 
itive charge. 
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Fig. 11. Optically thin gyrosynchrotron spectrum emitted in z 
direction. Different curves represent different times. O = 1 cor- 
responds to about 1 GHz. 

diet the emission frequency of the plasma waves. However, we 
may consult observed density profiles at height H ~ 3 10'* km, 
yielding n^ ~ 2.5 ■ 10^ cm"-^ averaged over the quiet corona 
(Fludra[l999 1. This number varies by a factor 3 when equator- 
to-pole variation is taken into account (Gallagher et al. 1999 1, 
and may be larger by a factor 10 above faculae (Dumont et 
al. [19821 Thus we associate a plasma density of 10^ cm"^ to 
2-10'' cm"^ with the z - H boundary of our simulation slab, 
corresponding to a plasma frequency of 100 to 400 MHz. 

4.4. Gyrosynchrotron emission 

The electrons which are trapped in closed loops are thought 
to be responsible for gyrosynchrotron radiation, manifesting in 
microwave solar continuum bursts'*. In this case there are many 
bursts which are optically thin at all frequencies (Fleishman et 
al. .2003j . and we assume here optically thin radiation for sim- 
pUcity. Moreover, we neglect the plasma response and consider 
emission in vacuum. Therefore, our results apply to the emis- 
sion at frequencies above the spectral peak provided by either 
optical thickness or the Razin effect. An electron on a circular 
orbit emits then at frequency nco/j the intensity (Schott. ri912> 



2 2 

n oj 



In = ^2~ (t^"^ ^ ■'ninv COS 0) + V^J'„^(m' COS 6)) 



r 



(15) 



where n is an integer, oj = eB/m is the local cyclotron fre- 
quency, and 6 is the angle between the gyration plane and the 
line of sight, which latter is taken along the z-direction. In our 
simulation, most trapped electrons have velocities v <k 1 (Fig. 
|6lbottom), so that the argument of the Bessel functions in Eq. 
J15> is small and J„{z) can be approximated by (2nnY^I^(ezl2)" 
(Abramowitz & Stegun. [1970> . yielding 



/«- 



w 



^ 1 + sin^e n |evcos0|2" 



(16) 



y^ cos^ 9 2n\ 2 i ' 

' In contrast, narrowband decimetric continua are better explained 
by transition radiation (Fleishman et al. l2005l Nita et al. l2b05 ). 



where Euler's number e - 2.718 is not to be confused with 
the elementary charge. The relative accuracy of the approxima- 
tion ( I16> is of order |vcos0| for all «; a numerical illustration 
is shown in Fig. ^| From Eq. ( I16> we see that the intensity 
/„ decays exponentially with n, so that the radiation is concen- 
trated at low harmonics. In the limits v — » and — » 90° 
only the fundamental {n - \) contributes. If the velocity has a 
component parallel to the magnetic field and to the observing 
direction, this results in a Doppler shift u) — > cf , and v in Eqns. 
( I15II16> has to be replaced by v^, the component perpendicu- 
lar to the magnetic field. Neglecting dispersion and absorption, 
we may thus obtain the radio spectrum by accumulating a his- 
togram of the frequencies nuf ly with weight /„ for all parti- 
cles and times. The result is shown in Figure^J Note that the 
spectrum decays monotonically and approximately exponen- 
tially with frequency, which is different from the familiar syn- 
chrotron (y » 1) shape with powerlaw rise and -decay. Here, 
the exponential decay is mostly caused by the harmonics de- 
pendence (Ea. ll6t but also supported by the exponential decay 
of magnetic field with height. Individual harmonics are washed 
out because of the magnetic field inhomogeneity. 
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N, /[km] Hjl al t? [c^fi] hajl 



1) 


5000 


128 


875 


40 


0.05 


4.38 


0.56 


2) 


5000 


256 


660 


50 


0.05 


6.84 


0.50 


3) 


5000 


256 


660 


50 


0.03 


6.84 


0.19 


4) 


5000 


256 


660 


64 


0.01 


6.84 


0.13 


5) 


1000 


256 


660 


64 


0.01 


6.84 


0.05 


6) 


1000 


256 


660 


64 


0.01 


1.37 


0.08 



Table 2. Simulation parameters. In all simulations, O ' = 
10"''s and uj = 0.2. 



4.5. Parameter exploration 

We have performed various production runs, involving a to- 
tal CPU time of about 4 months. Different runs have different 
values of the parameters H (slab height), a (force-free param- 
eter), and 77; see Table|2 The critical twist is fixed at uj - 0.2. 
The simulation discussed in Sections [4. II - 14. 4l is labeled 3 and 
represents an intermediate case. Each simulation involves sev- 
eral 10^ particles. An overview on the simulation results is pre- 
sented in Fig.^l The column a) contains the electron exit rates 
at z = (black line) and z-H (gray line) during the first third 
of the simulation time, normalized by the total number of sim- 
ulated electrons. The column b) shows the energy distributions 
at the end of the simulation (f - exit time, or t - r„iax)- The 
translation of the electron results onto HXR and radio wave 
predictions proceeds similar as in Sections [4. 2l to l4~4l and pre- 
serves the hardness ordering of the spectra. 

The simulations 1 to 6 may be summarized by saying that 
the results of Sections |4. 1 1 to l4~4l are generic, in the sense that 
the general shape and relative timing of the light curves are 
robust against variation of (H, a, 77). However, the time scal- 
ing, the division into precipitating and escaping populations, 
and the spectrum depend on the parameters (Tab.|2}. The time 
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Fig. 12. Parameter exploration. Left column: exit rates. Right column: terminal energy spectra. Black (gray, light gray) lines 
refers to precipitating (escaping, trapped) populations. See Tab.|2lfor the simulation parameters. 



delay between the precipitation pulse and the onset of escape 
(FigElleft) increases, naturally, with the slab height H. It also 
depends on the force-free parameter a, yet in a more subtle 
way. From the comparison of simulations 2 and 3 it is seen that 
a larger a favours prompt escape. This is not immediately ob- 
vious, since a larger a implies a smaller force-free scale 1 1 a 
and thus smaller magnetic loops which do not protrude to z = 
H. However, a larger a also admits more non-decaying modes 
(^j^ < a) and thus more open field lines along which the elec- 
trons can escape. Moreover, a larger a implies larger curvature, 
so that the criticality condition (Eq. |6j is more frequently met 
and the accelerating volume is larger {hd - 0.5/ in simulation 2 
but hd - 0.19/ in simulation 3). As a net effect, a large force- 
free parameter a thus favours rapid acceleration and escape. A 
similar enhancement of acceleration by the presence of small 
scales accounts for the difference between panels 4 and 5 in 
Fig -El These simulations differ only by the number of Fourier 



components. Run 4 contains wave vectors up to kj^l < 1.94, 
whereas run 5 contains wavevectors kj^l < 0.8. Accordingly, 
run 4 reaches somewhat higher energies (panels 4b and 5b). 
The escape time also increases (and the escape probability de- 
creases) with decreasing anomalous resistivity and thus with 
decreasing electric field (5a vs. 6a). The trapped component 
(light gray) is throughout softest, except for the case of very 
small resistivity (run 6) where all 3 populations behave simi- 
larly for E/mc^ S 0. 1 . Run 6 is similar to run 5 but with smaller 
resistivity, compensated by a somewhat larger dissipation vol- 
ume. 

As a general trend, it can be seen from Fig. E] (top to bot- 
tom) that smaller dissipation volumes and weaker electric fields 
yield a more gradual evolution and less intermediate energies 
(0.1 < E/mc^ < 1) but a comparable amount of high energies, 
so that the energy histograms have more tenuous but harder 
tails. This reflects a change of the nature of acceleration from 
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Fig. 13. The distribution of precipitating, escaping, plus trapped particle energy (black line) and the distribution of voltage drops 
(gray line). 



frequent and small energy increments to rare but violent energy 
gains, in the course of which the stochastic process leaves the 
domain of attraction of the central limit theorem. 

5. Summary and Discussion 

We have simulated gyrokinetic electron orbits in constant-a 
force-free magnetic fields with anomalous resistivity rj. The 
latter is localized in (postulated) dissipation regions where the 
magnetic twist V x b exceeds a given threshold Uc, correspond- 
ing to a critical scale l/uc ~ 3000 km. The dissipation regions 
cover about 10"^ of the simulated slab volume, and can be char- 
acterized by a column height h^ of 30 to 300 km. The (paral- 
lel) dissipative electric field exceeds the Dreicer field by one or 
two orders of magnitude and yields direct acceleration of run- 
away electrons. In general, the electrons visit numerous dissi- 
pation regions before reaching relativistic energies, and arrive 
at the chromosphere before escaping to the higher corona. The 
latter ordering agrees with an observed trend for HXR to pre- 



cede type III bursts (Aschwanden et al. 119921 Arzner & Benz 
12005 at . also in terms of absolute time delays. 

From a physics point of view, one may compare the en- 
ergies reached in Fig. E] to the voltage drop along the mag- 
netic field line inside a (simply connected) dissipation region 
D. (We use here the term 'voltage drop' rather than 'poten- 
tial drop' since E - rjj is not a potential field.) The quantity 
Ay = J E-dl gives an upper limit on the energy which can be 
gained inside a single dissipation region; reflection from con- 
verging field lines (and also deceleration from E itself) does 



generally prevent the particles from exploiting the full voltage 
drop. By choosing 5- 10^ random points inside the dissipation 
regions, and following the magnetic field lines going through 
these points, we find the distribution of voltage drops shown in 
FigEl(gr^y line), together with the distribution of all terminal 
electron energies (black line). Both kinetic energy and voltage 
drops are measured in units of the electron rest mass. As can 
be seen, the two distributions scale similarly in energy from 
one simulation to another (note the different energy axes!), and 
the sub-exponential tails of the particle energy occur above the 
largest voltage drop present in the simulation. This agrees with 
the observation that many dissipation regions are visited before 
the particles reach the highest energies. The root mean square 
voltage drops AVmis are of the order of the voltage drop across 
the magnetogram resolution, rjal. The distribution of voltage 
drops decays slower than the distribution of particle energies, 
which relates to the fact that AVrms is merely an upper bound on 
the kinetic energy gain. The peak of the particle distributions at 
lowest energies is mostly due to the trapped component. 

The present model predicts a statistical preference for the 
HXR to occur at one magnetic polarity, because of the assump- 
tion of constant a, by which the electric field E = rjaB is 
always parallel {a > 0) or antiparallel (a < 0) to the mag- 
netic field. The electrons preferentially impact the chromo- 
sphere not in the center of sunspots or filaments (where |B| is 
largest) but at their borders. Ions, due to their opposite charge, 
should preferentially impact at opposed magnetic polarity (but 
the ion dynamics is not directly comparable to the electron dy- 
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namics because of the large mass ratio). The experimental ver- 
ification or falsification of this model prediction, using non- 
randomized force-free extrapolations and spatially resolved X- 
ray (e.g., RHESSl) observations of flaring active regions, will 
be the subject of future investigations. It should also be pointed 
out that the impact map (Fig.|9j is obtained from simultaneous 
injection in all dissipation regions. If injection was restricted to 
a single dissipation region or group of dissipation regions, then 
only few impact regions would occur, as usually observed in 
individual solar flares. 

Finally, we should mention that the present approach has 
many caveats. Most prominent among them is the fact that test 
particle simulations can not tell us anything about the absolute 
number of electrons which are accelerated, and their applica- 
tion to the real solar corona requires caution in order to match 
global electrodynamic constraints (i.e. return currents; Spicer 
& Sudan 1984Jl. As a rule, our simulations can only account for 
a tenuous high-energy tail, because the back-action on the elec- 
tromagnetic field and Coulomb collisions are neglected. Then, 
we have used here magnetograms of non-flaring active regions 
- where the force-free extrapolation should be a good approx- 
imation - and have thus envisaged 'quiescent' coronal heating 
by nanoflares rather than large isolated events. The observable 
predictions, though, address sizable flares where all types of 
emission can be detected. We suggest that the actual flares are 
generated by a similar process, and should thus have similar 
characteristics. Conversely, the non-flaring active regions are 
constantly doing what flaring active regions do but in smaller, 
possibly undetected numbers. This coronal process may ac- 
count for spatial X-ray and radio fine structures in quiet solar 
regions (Benz et al. ll997> and non-thermal electrons measured 
in space during quiet times ('superhalo'; Lin 1998), as an al- 
ternative explanation of MHD wave acceleration in the solar 
wind. Also, the assumption of a time-independent magnetic 
field is an approximation only. It is certainly violated during 
large flares with re-structuring of the global magnetic topology. 
Also, it does not formally account for the continuous footpoint 
motion driving the Parker (1983 1 scenario. However, it may 
be a reasonable approximation for medium-size flares, where 
the Parker mechanism does not destroy the overall magnetic 
field structure. This point of view is supported by the work of 
Aulanier et al. (I2005all2005 b) who numerically study the evo- 
lution of magnetic flux tubes and quadrupolar configurations 
under sub-Alfvenic photospheric motion. Their (incompress- 
ible, resistive) MHD simulations demonstrate that current sheet 
formation and topological changes may be described by a se- 
quence of quasi-equilibria, so that a static approximation over 
a few ten seconds seems appropriate. With regard to particle 
acceleration, the temporal evolution of the magnetic fields is 
expected to decrease trapping and facilitate both precipitation 
and escape. 
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Fig. 14. Benchmark of the orbit integration schemes, showing 
position (left column) and velocity (right column) of a relativis- 
tic electron in strongly curved (///o=10) magnetic fields. SoUd 
line represents exact orbits (Eqns. l8l9t : dotted line, gyrokinetic 
approximation (Eqns.[TUI-ll4>. Gray crosses represent the hy- 
brid mode, with the exact regime marked dark and the gyroki- 
netic regime marked light. 



Appendix A: Numerical Implementation 

A.1. Regime switching and accuracy control 

The choice of the orbit integration method, exact or gyroki- 
netic, depends on two constraints. From a physical point of 
view, the gyrokinetic approximation requires the Larmor ra- 
dius to be small compared to the scale of the magnetic field, 
p <K /. From a numerical point of view, the gyrokinetic approx- 
imation is faster than exact orbit integration whenever the latter 
is dominated by particle motion rather than by field evaluation. 
Indeed, if t/ denotes the CPU time needed for a single field 
evaluation and Tg is the CPU time needed to resolve a single 
gyration, then the gyrokinetic approximation is by a factor 



T„ I 

Tf p 



(A.l) 



faster than exact orbit tracking. The time Tg is approxi- 
mately independent of the integration scheme as long as this 
is optimal, and a typical value on a medium-size workstation is 
found to be T^ ~ 5 ■ 10"^ s. The time cost of a single field eval- 
uation, on the same hardware, is tj ~ N^ x 5 ■ lO^^s with A^^ 
the number of Fourier components. Therefore the gyrokinetic 
method is numerically beneficial for A^^^ < 10 X I /p. Thus, if 
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the gyrokinetic approximation is physically allowed {lip » 1) 
then it is also numerically beneficial in the sparse-Fourier rep- 
resentation; the upper bound on A^^ is in practice rarely met. 

We turn now to the precise (and implemented) formula- 
tion of regime switching and accuracy control. The gyrokinetic 
approximation requires that the magnetic field should change 
slowly across a Larmor radius, (VB) ■ p < egB, where VB is the 
Jacobian and e^ <K 1 . For practical purposes, we replace this 
by the stronger and computationally less expensive constraint 
P^H/jl^y-S)!^ < e^B^ with p - \p\/B. The quantity p also in- 
cludes parallel momentum, so that the condition accounts for 
particle acceleration inside the dissipation regions. In addition 
to the gyrokinetic approximation error, there are numerical er- 
rors from the finite time step and field computation rate. These 
are controlled by limiting the Cash-Karp position- and momen- 
tum errors, and by enforcing re-calculation of the electromag- 
netic fields if |Axp 2/; \diBj\^ > e^, where Ax = x - Xom is 
the distance to the last field evaluation point Xo^. Using the 
same Jacobian, the fields are linearly extrapolated fromx^w. In 
order to avoid relentless switching between exact and gyroki- 
netic regimes, a minimum duration of one gyro period in each 
of the regimes is enforced. 



A.2. Benchmarks 

The exact orbits have been computed with diff^erent integration- 
and field evaluation schemes. The time integrators include the 
Verlet scheme and other leapfrog variants, and Runge-Kutta 
schemes with and without adaptive time stepping. The field 
evaluation was done either pointwise or with linear interpola- 
tion involving VB. Once the exact orbits were established, they 
have been used to benchmark the gyrokinetic orbits. As an ex- 
ample. Figure [2] displays individual cartesian components of 
position and velocity in an extreme situation, i.e. for a highly 
relativistic (v - 0.99c) particle with small magnetic scales (l/lo 
= 10). Solid line represents exact orbits, dotted line represents 
the gyrokinetic approximation, and gray represents the hybrid 
mode with automatic regime switching, where exact (gyroki- 
netic) regimes are marked dark (light). In the hybrid method, 
the particle starts with exact orbit integration and switches to 
the gyrokinetic description at Q.t ~ 10. It remains then in the 
gyrokinetic mode until at Qf ~ 120 a region of strong mag- 
netic curvature is encountered, where exact orbit tracing is en- 
forced (dark gray). At Qf ~ 200 the orbit switches back to the 
gyrokinetic method (light gray). During the gyrokinetic phase, 
the magnetic moment is conserved to within 0.1%, but varies 
during the exact-orbit phase by some 80%. As can be seen, the 
hybrid orbit is closer to the exact result than the purely gyroki- 
netic orbit. 
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